function r = hydro(Ti, TIME_STEPS)
  close all;
  fig = figure; 
  load('cm', 'mycmap');
  set(fig, 'Colormap', mycmap);


  GRID_SIZE = 128;
  DX = 0.01;
  DT = 0.0005;
  [gridx, gridy] = meshgrid([1:GRID_SIZE], [1:GRID_SIZE]);

  VIDEO = true;
  FPS = 30;
  SECONDS_PER_LIGHTYEAR = 10;
  if (VIDEO)
    vt = VideoWriter('v_t.avi', 'Archival');
    vvx = VideoWriter('v_vx.avi', 'Archival');
    vvy = VideoWriter('v_vy.avi', 'Archival');
    vv = VideoWriter('v_vv.avi', 'Archival');
    open(vt);
    open(vvx);
    open(vvy);
    open(vv);
  end

  ZERO_HYDRO = true;


  A = sqrt(3)*pi/18 - log(3)/2;

  vx = zeros(GRID_SIZE, GRID_SIZE);
  vy = zeros(GRID_SIZE, GRID_SIZE);
  px = zeros(GRID_SIZE, GRID_SIZE);
  py = zeros(GRID_SIZE, GRID_SIZE);
  Pi = zeros(GRID_SIZE, GRID_SIZE);
  B  = zeros(GRID_SIZE, GRID_SIZE);
  T  = ones(GRID_SIZE, GRID_SIZE);
  % shortcuts
  wx = zeros(GRID_SIZE, GRID_SIZE); % 1 - vx^2
  wy = zeros(GRID_SIZE, GRID_SIZE); % 1 - vy^2
  G = zeros(GRID_SIZE, GRID_SIZE); % sqrt(-1 + wx + wy)
  
  dtT = zeros(GRID_SIZE, GRID_SIZE); 
  dxdtT = zeros(GRID_SIZE, GRID_SIZE); 
  dydtT = zeros(GRID_SIZE, GRID_SIZE); 
  dxdxT = zeros(GRID_SIZE, GRID_SIZE); 
  dxdyT = zeros(GRID_SIZE, GRID_SIZE); 
  dydyT = zeros(GRID_SIZE, GRID_SIZE); 
  dtdtvx = zeros(GRID_SIZE, GRID_SIZE); 
  dtdtvy = zeros(GRID_SIZE, GRID_SIZE); 

  % derviatives
  dxvx = zeros(GRID_SIZE, GRID_SIZE);
  dxvy = zeros(GRID_SIZE, GRID_SIZE);
  dxpx = zeros(GRID_SIZE, GRID_SIZE);
  dxpy = zeros(GRID_SIZE, GRID_SIZE);
  dxPi = zeros(GRID_SIZE, GRID_SIZE);
  dxB  = zeros(GRID_SIZE, GRID_SIZE);
  dxT  = zeros(GRID_SIZE, GRID_SIZE);
  dyvx = zeros(GRID_SIZE, GRID_SIZE);
  dyvy = zeros(GRID_SIZE, GRID_SIZE);
  dypx = zeros(GRID_SIZE, GRID_SIZE);
  dypy = zeros(GRID_SIZE, GRID_SIZE);
  dyPi = zeros(GRID_SIZE, GRID_SIZE);
  dyB  = zeros(GRID_SIZE, GRID_SIZE);
  dyT  = zeros(GRID_SIZE, GRID_SIZE);
  dxdxvx = zeros(GRID_SIZE, GRID_SIZE);
  dxdyvx = zeros(GRID_SIZE, GRID_SIZE);
  dxdxvy = zeros(GRID_SIZE, GRID_SIZE);
  dxdyvy = zeros(GRID_SIZE, GRID_SIZE);
  dydyvx = zeros(GRID_SIZE, GRID_SIZE);
  dydyvy = zeros(GRID_SIZE, GRID_SIZE);


  % ********************************************
  % parallel computing
  % ********************************************



  sumT = 0;
  % initial conditions for vx vy come here
  % vx
  % vy
  % T
  TEMP_OFFSET = 100;
 T(3,3) = 0.8;
 cx = 0.5 * GRID_SIZE + 0.5;
 cy = 0.5 * GRID_SIZE + 0.5;
 for x = 1:GRID_SIZE
   for y = 1:GRID_SIZE
     T(x, y) = 0.4 * exp(-((x-cx)^2 + (y-cy)^2)/72);
   end
 end
  %T(:,:) = Ti(:,:);
  T = T + TEMP_OFFSET;
  for x = 1:GRID_SIZE
     for y = 1:GRID_SIZE
       % T(x,y) = T(x,y) +  0.01 * rand();
       %vx(x,y) = 0.001 * rand();
       %vy(x,y) = 0.001 * rand();
     end
  end

  % vx = 0.5 * ones(GRID_SIZE, GRID_SIZE);;
  vy = 0.9999999999* ones(GRID_SIZE, GRID_SIZE);;
  %****************************************************************
  % hydrolize the initial conditions
  %****************************************************************
  if (false)
  k = zeros(GRID_SIZE+1, GRID_SIZE+1);  
  [cx, cy] = meshgrid(0:GRID_SIZE-1,0:GRID_SIZE-1);
  [kx, ky] = meshgrid(-GRID_SIZE/2:GRID_SIZE/2);
  for kx = -GRID_SIZE/32:GRID_SIZE/32
    for ky = -GRID_SIZE/32:GRID_SIZE/32
      ft = T .* exp(-2 * pi * i * ( kx * cx + ky * cy ) / GRID_SIZE);
      k(kx+GRID_SIZE/2+1, ky+GRID_SIZE/2+1) = sum(sum(ft));
    end
  end
  T = zeros(GRID_SIZE:GRID_SIZE);
  for kx = -GRID_SIZE/32:GRID_SIZE/32
    for ky = -GRID_SIZE/32:GRID_SIZE/32
      T = T + k(kx+GRID_SIZE/2+1,ky+GRID_SIZE/2+1) * exp(2 * pi * i * ( kx * cx + ky * cy ) / GRID_SIZE);
    end
  end
  T = real(T / GRID_SIZE / GRID_SIZE);
end



  % ****************************************************************
  %  compute derivative kernels
  % ****************************************************************

  d1 = zeros(GRID_SIZE, GRID_SIZE);
  d2 = zeros(GRID_SIZE, GRID_SIZE);
  for x = 1:GRID_SIZE
    for y=1:GRID_SIZE
      s1 = 0;
      s2 = 0;
      for k = -GRID_SIZE/2 : GRID_SIZE/2
        s1 = s1 + 2. * pi * i / GRID_SIZE / GRID_SIZE * k * exp(2 * pi * i * k * (x - y) / GRID_SIZE);
        s2 = s2 + (2. * pi * i / GRID_SIZE * k)^2 / GRID_SIZE * exp(2 * pi * i * k * (x - y) / GRID_SIZE);
      end
      d1(x,y) = real(s1);
      d2(x,y) = real(s2);
    end
  end

  %dxT = dxdx(T);
  %dxT2 = dxdx2(T);
  % image(imagesc(dxT));
  %figure;
  %image(imagesc(dxT2));
  %image((dxT) * 128  + 128); 
  % image((dxT) * 128 + 128); 
  % truesize(fig);
  %fig=figure;
  %truesize(fig);
  %keyboard();
 % image((T-TEMP_OFFSET)*256+128);

  %return

  % ****************************************************************
  % compute the higher order initial conditions from ideal hydro:
  % ****************************************************************

  ts = 1;
  if (ZERO_HYDRO)
    ts = TIME_STEPS;
  end;

  tic;
  for t=1:ts

  wx = 1 - vx .* vx;
  wy = 1 - vy .* vy;
  G = sqrt(-1 + wx + wy);


  % compute the gradients
  dxvx = dx(vx);
  dyvx = dy(vx);
  dxvy = dx(vy);
  dyvy = dy(vy);
  dxT = dx(T);
  dyT = dy(T);
  if (~ZERO_HYDRO)
    dxdxvx = dxdx(vx);
    dxdyvx = dx(dy(vx));
    dydyvx = dydy(vx);

    dxdxvy = dxdx(vy);
    dxdyvy = dx(dy(vy));
    dydyvy = dydy(vy);

    dxdxT = dxdx(T);
    dxdyT = dx(dy(T));
    dydyT = dydy(T);
  end


      px = -(1./((wx+wy).*T)).*((wx+wy) .* T .* vy .* dyvx+(2 .* wx.^2+3.*wx .* (-1+wy)+(-1+wy).^2) .* dxT-vx.*((-1+wx+wy) .* vy .* dyT+T .* ((-1+wx+wy) .* dyvy-dxvx)));
      py = -((1./((wx + wy).*T)).*((1 + wx.^2 - 3.*wy + 2.*wy.^2 + wx.*(-2 + 3.*wy)).*dyT - (-1 + wx + wy).*vx.*vy.*dxT + T.*(vy.*(dyvy - (-1 + wx + wy).*dxvx) + (wx + wy).*vx.*dxvy)));
      dtT = -((vy.*dyT + vx.*dxT + T.*(dyvy+ dxvx))./(wx + wy));

  if (false)
  dxpx = dx(px);
  dypx = dy(px);
  dxpy = dx(py);
  dypy = dy(py);
  dxdtT = dx(dtT);
  dydtT = dy(dtT);

  for x = 1:GRID_SIZE
    for y = 1:GRID_SIZE
      % dtdtvx(x, y) = (1/((wx(x,y) + wy(x,y))^2*T(x,y)^2))*      (dxT(x,y)*((-1 + wx(x,y) + wy(x,y))*(wx(x,y) + wy(x,y))*                (-1 + 2*wx(x,y) + wy(x,y))*dtT(x,y) +                     2*(-1 + wx(x,y) + wx(x,y)^2 + 2*wx(x,y)*wy(x,y) + wy(x,y)^2)*T(x,y)*vy(x,y)*                py(x,y)) -          T(x,y)*(T(x,y)*((-((wx(x,y) + wx(x,y)^2 + 2*wx(x,y)*wy(x,y) + (-2 + wy(x,y))*(1 + wy(x,y)))*                                            dyvy(x,y) + (-2 + wx(x,y) - wy(x,y))*                             dxvx(x,y)))*          px(x,y) + (wx(x,y) + wy(x,y))^2*          dyvx(x,y)*                     py(x,y)) + vy(x,y)*                ((-(wx(x,y) + wx(x,y)^2 + 2*wx(x,y)*wy(x,y) + (-2 + wy(x,y))*(1 + wy(x,y))))*                     dyT(x,y)*          px(x,y) +                    (wx(x,y) + wy(x,y))^2*T(x,y)*          dypx(x,y)) +               (-1 + wx(x,y) + wy(x,y))*(wx(x,y) + wy(x,y))*(-1 + 2*wx(x,y) + wy(x,y))*       dxdtT(x,y)) + vx(x,y)*((-vy(x,y))*                (2*T(x,y)^2*(dyvy(x,y) +                         dxvx(x,y))*          py(x,y) + (-1 + wx(x,y) + wy(x,y))*(wx(x,y) +             wy(x,y))*(dyT(x,y)*                          dtT(x,y) - T(x,y)*                          dydtT(x,y))) +       T(x,y)*                (2*(-1 - wy(x,y) + 2*(wx(x,y) + wy(x,y))^2)*          dxT(x,y)*                               px(x,y) + (-2 + (-1 + wx(x,y))*wx(x,y) + wy(x,y) +                         2*wx(x,y)*wy(x,y) + wy(x,y)^2)*          dyT(x,y)*                     py(x,y) + (wx(x,y) + wy(x,y))*          T(x,y)*                     ((-1 + wx(x,y) + wy(x,y))*             dypy(x,y) -                        dxpx(x,y)))));

      % dtdtvy(x,y) = (1/((wx(x,y) + wy(x,y))^2*T(x,y)^2))*      (-2*T(x,y)*vx(x,y)^2*vy(x,y)*    dxT(x,y)*           px(x,y) +    dyT(x,y)*           ((-1 + wx(x,y) + wy(x,y))*(wx(x,y) + wy(x,y))*(-1 + wx(x,y) + 2*wy(x,y))*       dtT(x,y) +       2*T(x,y)*((-1 + wy(x,y) + (wx(x,y) + wy(x,y))^2)*vx(x,y)*                               px(x,y) + (-1 + wx(x,y)*(-1 + 2*wx(x,y)) +                         4*wx(x,y)*wy(x,y) + 2*wy(x,y)^2)*vy(x,y)*          py(x,y))) +    T(x,y)*(-2*T(x,y)*vy(x,y)^2*                (dyvy(x,y) +          dxvx(x,y))*py(x,y) +               (wx(x,y) +          wy(x,y))*(T(x,y)*((-(wx(x,y) + wy(x,y)))*dxvy(x,y)*                          px(x,y) -                         (dyvy(x,y) - (-1 + wx(x,y) + wy(x,y))*                               dxvx(x,y))*             py(x,y)) - (-1 + wx(x,y) + wy(x,y))*(-1 + wx(x,y) + 2*wy(x,y))*          dydtT(x,y)) + (wx(x,y) + wy(x,y))*vy(x,y)*                ((-1 + wx(x,y) + wy(x,y))*dxT(x,y)*                     px(x,y) + T(x,y)*                     (-dypy(x,y) + (-1 + wx(x,y) + wy(x,y))*                          dxpx(x,y)))) -          vx(x,y)*(2*T(x,y)*vy(x,y)^2*dxT(x,y)*                py(x,y) + vy(x,y)*                (2*T(x,y)^2*(dyvy(x,y) +                         dxvx(x,y))*          px(x,y) + (-1 + wx(x,y) + wy(x,y))*(wx(x,y) +             wy(x,y))*(dxT(x,y)*                          dtT(x,y) - T(x,y)*                          dxdtT(x,y))) + (wx(x,y) +          wy(x,y))*T(x,y)*                ((-(-1 + wx(x,y) + wy(x,y)))*dxT(x,y)*                     py(x,y) + (wx(x,y) + wy(x,y))*          T(x,y)*                     dxpy(x,y))));



      dtdtvx(x,y) = (1/(T(x,y)^2*(-2 + vx(x,y)^2 + vy(x,y)^2)^3))*(vx(x,y)^9*(T(x,y)*dydyT(x,y) + 8*dxT(x,y)^2) - vx(x,y)^8*(-10*vy(x,y)*dyT(x,y)* dxT(x,y) + T(x,y)*(-7*dyvy(x,y)* dxT(x,y) + dyT(x,y)*(-dyvx(x,y) + dxvy(x,y)) + vy(x,y)*dxdyT(x,y)) + T(x,y)^2*dxdyvy(x,y)) + vx(x,y)^5*(-5*T(x,y)^2*dyvy(x,y)^2 + 13*T(x,y)*dydyT(x,y) + 112*dxT(x,y)^2 + 6*vy(x,y)^4*(dyT(x,y)^2 + T(x,y)*dydyT(x,y) + 7*dxT(x,y)^2) - 7*T(x,y)^2*dyvy(x,y)* dxvx(x,y) - 4*T(x,y)^2*dxvx(x,y)^2 - 2*T(x,y)^2*dyvx(x,y)* dxvy(x,y) + T(x,y)* vy(x,y)*((-dyT(x,y))*(25* dyvy(x,y) + 9*dxvx(x,y)) + 2*dxT(x,y)*(23*dyvx(x,y) + 6*dxvy(x,y)) + T(x,y)*(7*dydyvy(x,y) - 5*dxdyvx(x,y))) - 3*T(x,y)* vy(x,y)^3*((-dyT(x,y))*(5* dyvy(x,y) + dxvx(x,y)) + 2*dxT(x,y)*(4*dyvx(x,y) + dxvy(x,y)) + T(x,y)*(dydyvy(x,y) - dxdyvx(x,y))) + vy(x,y)^2*(-14*dyT(x,y)^2 - 142*dxT(x,y)^2 + 3*T(x,y)^2* dyvy(x,y)*(dyvy(x,y) + dxvx(x,y)) - 2*T(x,y)*(8*dydyT(x,y) - 5*dxdxT(x,y))) - 16*T(x,y)*dxdxT(x,y)) + vx(x,y)^7*(T(x,y)^2*dyvy(x,y)^2 - 6*T(x,y)*dydyT(x,y) - 52*dxT(x,y)^2 + 2*vy(x,y)^2*(dyT(x,y)^2 + 2*T(x,y)*dydyT(x,y) + 15*dxT(x,y)^2) + T(x,y)^2*dyvy(x,y)* dxvx(x,y) + T(x,y)* vy(x,y)*(dyT(x,y)*(5*dyvy(x,y) + dxvx(x,y)) - 2*dxT(x,y)*(4*dyvx(x,y) + dxvy(x,y)) + T(x,y)*(-dydyvy(x,y) + dxdyvx(x,y))) + 4*T(x,y)*dxdxT(x,y)) + vx(x,y)*(4*T(x,y)*dydyT(x,y) + 32*dxT(x,y)^2 + vy(x,y)^8*(2*dyT(x,y)^2 + T(x,y)*dydyT(x,y) + 6*dxT(x,y)^2) - 8*T(x,y)^2*dxvx(x,y)^2 - 8*T(x,y)^2*dyvx(x,y)* dxvy(x,y) + T(x,y)* vy(x,y)^5*((-dyT(x,y))*(25* dyvy(x,y) + 9*dxvx(x,y)) + 2*dxT(x,y)*(23*dyvx(x,y) + 6*dxvy(x,y)) + T(x,y)*(7*dydyvy(x,y) - 5*dxdyvx(x,y))) - 2*T(x,y)* vy(x,y)^3*(-2* dyT(x,y)*(7*dyvy(x,y) + 3*dxvx(x,y)) + dxT(x,y)*(43*dyvx(x,y) + 13*dxvy(x,y)) + T(x,y)*(7*dydyvy(x,y) - 5*dxdyvx(x,y))) + 4*T(x,y)* vy(x,y)*((-dyT(x,y))*(3* dyvy(x,y) + dxvx(x,y)) + dxT(x,y)*(13*dyvx(x,y) + 5*dxvy(x,y)) + 2*T(x,y)*(dydyvy(x,y) - dxdyvx(x,y))) + T(x,y)* vy(x,y)^7*(dyT(x,y)*(5*dyvy(x,y) + dxvx(x,y)) - 2*dxT(x,y)*(4*dyvx(x,y) + dxvy(x,y)) + T(x,y)*(-dydyvy(x,y) + dxdyvx(x,y))) + vy(x,y)^4*(24*dyT(x,y)^2 + 88*dxT(x,y)^2 - T(x,y)^2*(5*dyvy(x,y)^2 + 7*dyvy(x,y)* dxvx(x,y) + 4*dxvx(x,y)^2 + 2*dyvx(x,y)* dxvy(x,y)) + T(x,y)*(7*dydyT(x,y) - 10*dxdxT(x,y))) - 2*vy(x,y)^2*(6*dyT(x,y)^2 + 44*dxT(x,y)^2 - T(x,y)^2*(dyvy(x,y)* dxvx(x,y) + 5*dxvx(x,y)^2 + 4*dyvx(x,y)* dxvy(x,y)) + 4*T(x,y)*(dydyT(x,y) - 2*dxdxT(x,y))) - 8*T(x,y)*dxdxT(x,y) + vy(x,y)^6*(-14*dyT(x,y)^2 - 38*dxT(x,y)^2 + T(x,y)^2* dyvy(x,y)*(dyvy(x,y) + dxvx(x,y)) + T(x,y)*(-4*dydyT(x,y) + 2*dxdxT(x,y)))) + vx(x,y)^3*(vy(x,y)^6*(6*dyT(x,y)^2 + 4*T(x,y)*dydyT(x,y) + 26*dxT(x,y)^2) - 2*T(x,y)* vy(x,y)*(-4*dyT(x,y)*(4*dyvy(x,y) + dxvx(x,y)) + 15*dxT(x,y)*(3*dyvx(x,y) + dxvy(x,y)) + T(x,y)*(7*dydyvy(x,y) - 5*dxdyvx(x,y))) + 2*T(x,y)* vy(x,y)^3*((-dyT(x,y))*(25* dyvy(x,y) + 9*dxvx(x,y)) + 2*dxT(x,y)*(23*dyvx(x,y) + 6*dxvy(x,y)) + T(x,y)*(7*dydyvy(x,y) - 5*dxdyvx(x,y))) - 3*T(x,y)* vy(x,y)^5*((-dyT(x,y))*(5* dyvy(x,y) + dxvx(x,y)) + 2*dxT(x,y)*(4*dyvx(x,y) + dxvy(x,y)) + T(x,y)*(dydyvy(x,y) - dxdyvx(x,y))) + 2*vy(x,y)^2*(12*dyT(x,y)^2 + 100*dxT(x,y)^2 - T(x,y)^2*(5*dyvy(x,y)^2 + 7*dyvy(x,y)* dxvx(x,y) + 4*dxvx(x,y)^2 + 2*dyvx(x,y)* dxvy(x,y)) + T(x,y)*(10*dydyT(x,y) - 13*dxdxT(x,y))) + vy(x,y)^4*(-28*dyT(x,y)^2 - 128*dxT(x,y)^2 + 3*T(x,y)^2* dyvy(x,y)*(dyvy(x,y) + dxvx(x,y)) + T(x,y)*(-14*dydyT(x,y) + 8*dxdxT(x,y))) + 2*(-50*dxT(x,y)^2 + T(x,y)^2*(2*dyvy(x,y)^2 + dyvy(x,y)* dxvx(x,y) + 3*dxvx(x,y)^2 + 4*dyvx(x,y)* dxvy(x,y)) + T(x,y)*(-6*dydyT(x,y) + 10*dxdxT(x,y)))) - vx(x,y)^4*(-2* vy(x,y)*(68 - 85*vy(x,y)^2 + 24*vy(x,y)^4)* dyT(x,y)* dxT(x,y) + T(x,y)*((-59 + 68*vy(x,y)^2 - 18*vy(x,y)^4)* dyvy(x,y)* dxT(x,y) + 32*dxT(x,y)* dxvx(x,y) - 6*vy(x,y)^4*dxT(x,y)* dxvx(x,y) + dyT(x,y)*((2 - 16*vy(x,y)^2 + 9*vy(x,y)^4)* dyvx(x,y) + (19 - 15*vy(x,y)^2 + 3*vy(x,y)^4)* dxvy(x,y)) + 35*vy(x,y)*dxdyT(x,y) - 32*vy(x,y)^3*dxdyT(x,y) + 6*vy(x,y)^5*dxdyT(x,y)) + T(x,y)^2*(3*vy(x,y)^3* dyvx(x,y)*(dyvy(x,y) + dxvx(x,y)) - vy(x,y)*(dyvy(x,y) + dxvx(x,y))*(3* dyvx(x,y) - 2*dxvy(x,y)) - 3*vy(x,y)^4*(dydyvx(x,y) - dxdyvy(x,y)) + 19*dxdyvy(x,y) + vy(x,y)^2*(6*dydyvx(x,y) - 17*dxdyvy(x,y) - 5*dxdxvx(x,y)) + 7*dxdxvx(x,y))) - vx(x,y)^2*(2*vy(x,y)*(56 - 112*vy(x,y)^2 + 71*vy(x,y)^4 - 14*vy(x,y)^6)*dyT(x,y)* dxT(x,y) + T(x,y)*(2*(16 - 18*vy(x,y)^2 + 11*vy(x,y)^4 - 2*vy(x,y)^6)*dyvy(x,y)* dxT(x,y) - 36*dxT(x,y)* dxvx(x,y) + 22*vy(x,y)^2*dxT(x,y)* dxvx(x,y) + 12*vy(x,y)^4*dxT(x,y)* dxvx(x,y) - 6*vy(x,y)^6*dxT(x,y)* dxvx(x,y) + dyT(x,y)*((-12 + 48*vy(x,y)^2 - 41*vy(x,y)^4 + 11*vy(x,y)^6)* dyvx(x,y) + (-16 + 22*vy(x,y)^2 - 6*vy(x,y)^4 + vy(x,y)^6)*dxvy(x,y)) - 36*vy(x,y)*dxdyT(x,y) + 58*vy(x,y)^3*dxdyT(x,y) - 28*vy(x,y)^5*dxdyT(x,y) + 4*vy(x,y)^7*dxdyT(x,y)) + T(x,y)^2*(3*vy(x,y)^5* dyvx(x,y)*(dyvy(x,y) + dxvx(x,y)) - 2*vy(x,y)^3*(dyvy(x,y) + dxvx(x,y))*(3* dyvx(x,y) - 2*dxvy(x,y)) + 2*vy(x,y)*(dyvy(x,y) + dxvx(x,y))*(2* dyvx(x,y) - dxvy(x,y)) + vy(x,y)^6*(-3*dydyvx(x,y) + dxdyvy(x,y)) + 2*vy(x,y)^4*(6*dydyvx(x,y) - 5*dxdyvy(x,y) - 2*dxdxvx(x,y)) - 12*vy(x,y)^2*(dydyvx(x,y) - 2*dxdyvy(x,y) - dxdxvx(x,y)) - 8*(2*dxdyvy(x,y) + dxdxvx(x,y)))) - vx(x,y)^6*(6*vy(x,y)*(11 - 6*vy(x,y)^2)* dyT(x,y)* dxT(x,y) + T(x,y)*(dyT(x,y)*((3 + vy(x,y)^2)* dyvx(x,y) + (-8 + 3*vy(x,y)^2)* dxvy(x,y)) - 2*((-19 + 10*vy(x,y)^2)* dyvy(x,y)* dxT(x,y) + (2 + vy(x,y)^2)* dxT(x,y)* dxvx(x,y) - 2*vy(x,y)*(-3 + vy(x,y)^2)* dxdyT(x,y))) + T(x,y)^2*(vy(x,y)* dyvx(x,y)*(dyvy(x,y) + dxvx(x,y)) - vy(x,y)^2*(dydyvx(x,y) - 3*dxdyvy(x,y)) - 2*(4*dxdyvy(x,y) + dxdxvx(x,y)))) - (-2 + vy(x,y)^2)*(2* vy(x,y)*(8 - 18*vy(x,y)^2 + 13*vy(x,y)^4 - 3*vy(x,y)^6)*dyT(x,y)* dxT(x,y) + T(x,y)*((-1 + vy(x,y)^2)* dyT(x,y)*(2*(2 - 5*vy(x,y)^2 + 2*vy(x,y)^4)* dyvx(x,y) + (2 + vy(x,y)^2)* dxvy(x,y)) + (-2 + vy(x,y)^2)*((-1 - 4*vy(x,y)^2 + vy(x,y)^4)* dyvy(x,y)* dxT(x,y) + (-1 + vy(x,y)^2)*(-2*(1 + vy(x,y)^2)* dxT(x,y)* dxvx(x,y) + vy(x,y)*(-3 + vy(x,y)^2)* dxdyT(x,y)))) - T(x,y)^2*(vy(x,y)^6*dydyvx(x,y) - vy(x,y)^5* dyvx(x,y)*(dyvy(x,y) + dxvx(x,y)) + vy(x,y)^3*(dyvy(x,y) + dxvx(x,y))*(dyvx(x,y) - 2*dxvy(x,y)) + 2*vy(x,y)*(dyvy(x,y) + dxvx(x,y))*(dyvx(x,y) + dxvy(x,y)) + 2*(dxdyvy(x,y) + dxdxvx(x,y)) + vy(x,y)^4*(-4*dydyvx(x,y) + dxdyvy(x,y) + dxdxvx(x,y)) + vy(x,y)^2*(4*dydyvx(x,y) - 3*(dxdyvy(x,y) + dxdxvx(x,y))))));


      dtdtvy(x,y) =  (1/(T(x,y)^2*(-2 + vx(x,y)^2 + vy(x,y)^2)^3))*(vy(x,y)^8*(vx(x,y)*(10*dyT(x,y)* dxT(x,y) - T(x,y)*dxdyT(x,y)) + T(x,y)*((-dyvx(x,y))* dxT(x,y) + 7*dyT(x,y)* dxvx(x,y) + dxT(x,y)* dxvy(x,y) - T(x,y)*dxdyvx(x,y))) + vy(x,y)^9*(8*dyT(x,y)^2 + T(x,y)*dxdxT(x,y)) + vy(x,y)^7*((-52 + 30*vx(x,y)^2)* dyT(x,y)^2 + 2*vx(x,y)^2*dxT(x,y)^2 - 2*T(x,y)*vx(x,y)* dyT(x,y)*(dyvx(x,y) + 4*dxvy(x,y)) + T(x,y)*(4*dydyT(x,y) + vx(x,y)* dxT(x,y)*(dyvy(x,y) + 5*dxvx(x,y)) - 6*dxdxT(x,y) + 4*vx(x,y)^2*dxdxT(x,y)) + T(x,y)^2*(dyvy(x,y)* dxvx(x,y) + dxvx(x,y)^2 + vx(x,y)*(dxdyvy(x,y) - dxdxvx(x,y)))) + vy(x,y)^5*(2*(56 - 71*vx(x,y)^2 + 21*vx(x,y)^4)* dyT(x,y)^2 + 2*vx(x,y)^2*(-7 + 3*vx(x,y)^2)* dxT(x,y)^2 + 2*T(x,y)*vx(x,y)* dyT(x,y)*(-3*(-2 + vx(x,y)^2)* dyvx(x,y) + (23 - 12*vx(x,y)^2)* dxvy(x,y)) + T(x,y)*(2*(-8 + 5*vx(x,y)^2)* dydyT(x,y) + 3*vx(x,y)^3* dxT(x,y)*(dyvy(x,y) + 5*dxvx(x,y)) - vx(x,y)* dxT(x,y)*(9*dyvy(x,y) + 25*dxvx(x,y)) + 13*dxdxT(x,y) - 16*vx(x,y)^2*dxdxT(x,y) + 6*vx(x,y)^4*dxdxT(x,y)) - T(x,y)^2*(4* dyvy(x,y)^2 + (7 - 3*vx(x,y)^2)* dyvy(x,y)* dxvx(x,y) + (5 - 3*vx(x,y)^2)* dxvx(x,y)^2 + 2*dyvx(x,y)* dxvy(x,y) + 5*vx(x,y)*dxdyvy(x,y) - 3*vx(x,y)^3*dxdyvy(x,y) - 7*vx(x,y)*dxdxvx(x,y) + 3*vx(x,y)^3*dxdxvx(x,y))) + vy(x,y)^3*(2*(-50 + 100*vx(x,y)^2 - 64*vx(x,y)^4 + 13*vx(x,y)^6)*dyT(x,y)^2 + 2*vx(x,y)^2*(12 - 14*vx(x,y)^2 + 3*vx(x,y)^4)* dxT(x,y)^2 - 2*T(x,y)*vx(x,y)* dyT(x,y)*(3*(5 - 4*vx(x,y)^2 + vx(x,y)^4)* dyvx(x,y) + (45 - 46*vx(x,y)^2 + 12*vx(x,y)^4)* dxvy(x,y)) + T(x,y)*((20 - 26*vx(x,y)^2 + 8*vx(x,y)^4)* dydyT(x,y) + 8*vx(x,y)* dxT(x,y)*(dyvy(x,y) + 4*dxvx(x,y)) + 3*vx(x,y)^5* dxT(x,y)*(dyvy(x,y) + 5*dxvx(x,y)) - 2*vx(x,y)^3* dxT(x,y)*(9*dyvy(x,y) + 25*dxvx(x,y)) - 12*dxdxT(x,y) + 20*vx(x,y)^2*dxdxT(x,y) - 14*vx(x,y)^4*dxdxT(x,y) + 4*vx(x,y)^6*dxdxT(x,y)) + T(x,y)^2*((6 - 8*vx(x,y)^2)* dyvy(x,y)^2 + (2 - 14*vx(x,y)^2 + 3*vx(x,y)^4)* dyvy(x,y)* dxvx(x,y) + (4 - 10*vx(x,y)^2 + 3*vx(x,y)^4)* dxvx(x,y)^2 + 8*dyvx(x,y)* dxvy(x,y) - 4*vx(x,y)^2*dyvx(x,y)* dxvy(x,y) + 10*vx(x,y)*dxdyvy(x,y) - 10*vx(x,y)^3*dxdyvy(x,y) + 3*vx(x,y)^5*dxdyvy(x,y) - 14*vx(x,y)*dxdxvx(x,y) + 14*vx(x,y)^3*dxdxvx(x,y) - 3*vx(x,y)^5*dxdxvx(x,y))) + vy(x,y)*(2*(-2 + vx(x,y)^2)^2*(4 - 7*vx(x,y)^2 + 3*vx(x,y)^4)*dyT(x,y)^2 + 2*vx(x,y)^2*(-6 + 12*vx(x,y)^2 - 7*vx(x,y)^4 + vx(x,y)^6)*dxT(x,y)^2 - 2*T(x,y)*vx(x,y)*(-2 + vx(x,y)^2)* dyT(x,y)*((5 - 4*vx(x,y)^2 + vx(x,y)^4)* dyvx(x,y) + (13 - 15*vx(x,y)^2 + 4*vx(x,y)^4)* dxvy(x,y)) + T(x,y)*(2*(-2 + vx(x,y)^2)^2*(-1 + vx(x,y)^2)* dydyT(x,y) - 4*vx(x,y)* dxT(x,y)*(dyvy(x,y) + 3*dxvx(x,y)) + vx(x,y)^7* dxT(x,y)*(dyvy(x,y) + 5*dxvx(x,y)) + 4*vx(x,y)^3* dxT(x,y)*(3*dyvy(x,y) + 7*dxvx(x,y)) - vx(x,y)^5* dxT(x,y)*(9*dyvy(x,y) + 25*dxvx(x,y)) + 4*dxdxT(x,y) - 8*vx(x,y)^2*dxdxT(x,y) + 7*vx(x,y)^4*dxdxT(x,y) - 4*vx(x,y)^6*dxdxT(x,y) + vx(x,y)^8*dxdxT(x,y)) + T(x,y)^2*(-2*(4 - 5*vx(x,y)^2 + 2*vx(x,y)^4)* dyvy(x,y)^2 + vx(x,y)^2*(2 - 7*vx(x,y)^2 + vx(x,y)^4)* dyvy(x,y)* dxvx(x,y) + vx(x,y)^6*dxvx(x,y)^2 - 8*dyvx(x,y)* dxvy(x,y) + 8*vx(x,y)^2*dyvx(x,y)* dxvy(x,y) - vx(x,y)^4*(5*dxvx(x,y)^2 + 2*dyvx(x,y)* dxvy(x,y)) + 2*vx(x,y)^3*(5*dxdyvy(x,y) - 7*dxdxvx(x,y)) - 8*vx(x,y)*(dxdyvy(x,y) - dxdxvx(x,y)) + vx(x,y)^7*(dxdyvy(x,y) - dxdxvx(x,y)) + vx(x,y)^5*(-5*dxdyvy(x,y) + 7*dxdxvx(x,y)))) - vy(x,y)^2*(2*vx(x,y)*(56 - 112*vx(x,y)^2 + 71*vx(x,y)^4 - 14*vx(x,y)^6)*dyT(x,y)* dxT(x,y) + T(x,y)*((-16 + 22*vx(x,y)^2 - 6*vx(x,y)^4 + vx(x,y)^6)*dyvx(x,y)* dxT(x,y) + dyT(x,y)*((-36 + 22*vx(x,y)^2 + 12*vx(x,y)^4 - 6*vx(x,y)^6)*dyvy(x,y) + 2*(16 - 18*vx(x,y)^2 + 11*vx(x,y)^4 - 2*vx(x,y)^6)*dxvx(x,y)) - 12*dxT(x,y)* dxvy(x,y) + 48*vx(x,y)^2*dxT(x,y)* dxvy(x,y) - 41*vx(x,y)^4*dxT(x,y)* dxvy(x,y) + 11*vx(x,y)^6*dxT(x,y)* dxvy(x,y) - 36*vx(x,y)*dxdyT(x,y) + 58*vx(x,y)^3*dxdyT(x,y) - 28*vx(x,y)^5*dxdyT(x,y) + 4*vx(x,y)^7*dxdyT(x,y)) + T(x,y)^2*(2* vx(x,y)^3*(dyvy(x,y) + dxvx(x,y))*(2* dyvx(x,y) - 3*dxvy(x,y)) - 2*vx(x,y)*(dyvy(x,y) + dxvx(x,y))*(dyvx(x,y) - 2*dxvy(x,y)) + 3*vx(x,y)^5*(dyvy(x,y) + dxvx(x,y))* dxvy(x,y) - 8*(dydyvy(x,y) + 2*dxdyvx(x,y)) - 2*vx(x,y)^4*(2*dydyvy(x,y) + 5*dxdyvx(x,y) - 6*dxdxvy(x,y)) + vx(x,y)^6*(dxdyvx(x,y) - 3*dxdxvy(x,y)) + 12*vx(x,y)^2*(dydyvy(x,y) + 2*dxdyvx(x,y) - dxdxvy(x,y)))) - vy(x,y)^6*(6*vx(x,y)*(11 - 6*vx(x,y)^2)* dyT(x,y)* dxT(x,y) + T(x,y)*((-8 + 3*vx(x,y)^2)* dyvx(x,y)* dxT(x,y) - 2*dyT(x,y)*((2 + vx(x,y)^2)* dyvy(x,y) + (-19 + 10*vx(x,y)^2)* dxvx(x,y)) + 3*dxT(x,y)* dxvy(x,y) + vx(x,y)^2*dxT(x,y)* dxvy(x,y) - 12*vx(x,y)*dxdyT(x,y) + 4*vx(x,y)^3*dxdyT(x,y)) + T(x,y)^2*(-2*dydyvy(x,y) + vx(x,y)*(dyvy(x,y) + dxvx(x,y))* dxvy(x,y) - 8*dxdyvx(x,y) + vx(x,y)^2*(3*dxdyvx(x,y) - dxdxvy(x,y)))) - (-2 + vx(x,y)^2)*(2* vx(x,y)*(8 - 18*vx(x,y)^2 + 13*vx(x,y)^4 - 3*vx(x,y)^6)*dyT(x,y)* dxT(x,y) + T(x,y)*((-(-2 + vx(x,y)^2))* dyT(x,y)*(2*(-1 + vx(x,y)^4)* dyvy(x,y) + (1 + 4*vx(x,y)^2 - vx(x,y)^4)* dxvx(x,y)) + (-1 + vx(x,y)^2)*((2 + vx(x,y)^2)* dyvx(x,y)* dxT(x,y) + (-2 + vx(x,y)^2)*(2*(-1 + 2*vx(x,y)^2)* dxT(x,y)* dxvy(x,y) + vx(x,y)*(-3 + vx(x,y)^2)* dxdyT(x,y)))) - T(x,y)^2*((-vx(x,y)^3)*(dyvy(x,y) + dxvx(x,y))*(2* dyvx(x,y) - dxvy(x,y)) - vx(x,y)^5*(dyvy(x,y) + dxvx(x,y))* dxvy(x,y) + 2*vx(x,y)*(dyvy(x,y) + dxvx(x,y))*(dyvx(x,y) + dxvy(x,y)) + 2*(dydyvy(x,y) + dxdyvx(x,y)) + vx(x,y)^4*(dydyvy(x,y) + dxdyvx(x,y) - 4*dxdxvy(x,y)) + vx(x,y)^6*dxdxvy(x,y) + vx(x,y)^2*(-3*dydyvy(x,y) - 3*dxdyvx(x,y) + 4*dxdxvy(x,y)))) - vy(x,y)^4*(-2* vx(x,y)*(68 - 85*vx(x,y)^2 + 24*vx(x,y)^4)* dyT(x,y)* dxT(x,y) + T(x,y)*((19 - 15*vx(x,y)^2 + 3*vx(x,y)^4)* dyvx(x,y)* dxT(x,y) + dyT(x,y)*((32 - 6*vx(x,y)^4)* dyvy(x,y) + (-59 + 68*vx(x,y)^2 - 18*vx(x,y)^4)* dxvx(x,y)) + 2*dxT(x,y)* dxvy(x,y) - 16*vx(x,y)^2*dxT(x,y)* dxvy(x,y) + 9*vx(x,y)^4*dxT(x,y)* dxvy(x,y) + 35*vx(x,y)*dxdyT(x,y) - 32*vx(x,y)^3*dxdyT(x,y) + 6*vx(x,y)^5*dxdyT(x,y)) + T(x,y)^2*(7*dydyvy(x,y) + vx(x,y)*(dyvy(x,y) + dxvx(x,y))*(2* dyvx(x,y) - 3*dxvy(x,y)) + 3*vx(x,y)^3*(dyvy(x,y) + dxvx(x,y))* dxvy(x,y) + 19*dxdyvx(x,y) + 3*vx(x,y)^4*(dxdyvx(x,y) - dxdxvy(x,y)) + vx(x,y)^2*(-5*dydyvy(x,y) - 17*dxdyvx(x,y) + 6*dxdxvy(x,y)))));



















% ********************************************


      Pi(x,y) = (1/(9*(-1 + wx(x,y) + wy(x,y))^(9/2)))*(pi*T(x,y)* (-8*pi*(-1 + wx(x,y) + wy(x,y))^2* T(x,y)*(wx(x,y)*wy(x,y)*dyvx(x,y) + wy(x,y)*(wx(x,y)*dxvy(x,y) + vy(x,y)*px(x,y)) + vx(x,y)*(vy(x,y)*(wx(x,y)*dyvy(x,y) + wy(x,y)*dxvx(x,y)) + wx(x,y)*py(x,y))) - 3*sqrt(-1 + wx(x,y) + wy(x,y))* (2*wx(x,y)*wy(x,y)*vy(x,y)*dydyvx(x,y) + 2*A*wx(x,y)*wy(x,y)*vy(x,y)*dydyvx(x,y) - 2*wx(x,y)^2*wy(x,y)*vy(x,y)*dydyvx(x,y) - 2*A*wx(x,y)^2*wy(x,y)*vy(x,y)*dydyvx(x,y) - 2*wx(x,y)*wy(x,y)^2*vy(x,y)*dydyvx(x,y) - 2*A*wx(x,y)*wy(x,y)^2*vy(x,y)*dydyvx(x,y) - 9*wx(x,y)*dyvy(x,y)* dxvy(x,y) - 11*A*wx(x,y)*dyvy(x,y)* dxvy(x,y) + 9*wx(x,y)^2*dyvy(x,y)*dxvy(x,y) + 11*A*wx(x,y)^2*dyvy(x,y)* dxvy(x,y) + 5*wx(x,y)*wy(x,y)*dyvy(x,y)* dxvy(x,y) + 7*A*wx(x,y)*wy(x,y)* dyvy(x,y)* dxvy(x,y) - 10*wx(x,y)^2*wy(x,y)*dyvy(x,y)* dxvy(x,y) - 10*A*wx(x,y)^2*wy(x,y)* dyvy(x,y)* dxvy(x,y) + 4*wx(x,y)*wy(x,y)^2*dyvy(x,y)* dxvy(x,y) + 4*A*wx(x,y)*wy(x,y)^2* dyvy(x,y)* dxvy(x,y) - 5*wy(x,y)*dxvx(x,y)* dxvy(x,y) - 3*A*wy(x,y)*dxvx(x,y)*dxvy(x,y) - 3*wx(x,y)*wy(x,y)*dxvx(x,y)* dxvy(x,y) - 5*A*wx(x,y)*wy(x,y)*dxvx(x,y)* dxvy(x,y) + 8*wx(x,y)^2*wy(x,y)* dxvx(x,y)* dxvy(x,y) + 8*A*wx(x,y)^2*wy(x,y)*dxvx(x,y)* dxvy(x,y) + 5*wy(x,y)^2*dxvx(x,y)*dxvy(x,y) + 3*A*wy(x,y)^2*dxvx(x,y)* dxvy(x,y) - 6*wx(x,y)*wy(x,y)^2*dxvx(x,y)* dxvy(x,y) - 6*A*wx(x,y)*wy(x,y)^2* dxvx(x,y)* dxvy(x,y) + 2*wx(x,y)*vy(x,y)*dxdyvy(x,y) + 2*A*wx(x,y)*vy(x,y)*dxdyvy(x,y) - 4*wx(x,y)^2*vy(x,y)*dxdyvy(x,y) - 4*A*wx(x,y)^2*vy(x,y)*dxdyvy(x,y) + 2*wx(x,y)^3*vy(x,y)*dxdyvy(x,y) + 2*A*wx(x,y)^3*vy(x,y)*dxdyvy(x,y) - 2*wx(x,y)*wy(x,y)^2*vy(x,y)*dxdyvy(x,y) - 2*A*wx(x,y)*wy(x,y)^2*vy(x,y)*dxdyvy(x,y) + 2*wy(x,y)*vy(x,y)*dxdxvx(x,y) + 2*A*wy(x,y)*vy(x,y)*dxdxvx(x,y) - 4*wx(x,y)*wy(x,y)*vy(x,y)*dxdxvx(x,y) - 4*A*wx(x,y)*wy(x,y)*vy(x,y)*dxdxvx(x,y) + 2*wx(x,y)^2*wy(x,y)*vy(x,y)*dxdxvx(x,y) + 2*A*wx(x,y)^2*wy(x,y)*vy(x,y)*dxdxvx(x,y) - 2*wy(x,y)^2*vy(x,y)*dxdxvx(x,y) - 2*A*wy(x,y)^2*vy(x,y)*dxdxvx(x,y) + 2*wx(x,y)*wy(x,y)^2*vy(x,y)*dxdxvx(x,y) + 2*A*wx(x,y)*wy(x,y)^2*vy(x,y)*dxdxvx(x,y) - 5*wx(x,y)*vy(x,y)*dyvy(x,y)* px(x,y) - 3*A*wx(x,y)*vy(x,y)* dyvy(x,y)* px(x,y) + 5*wx(x,y)^2*vy(x,y)*dyvy(x,y)* px(x,y) + 3*A*wx(x,y)^2*vy(x,y)* dyvy(x,y)* px(x,y) - 2*wy(x,y)*vy(x,y)*dyvy(x,y)* px(x,y) - 2*A*wy(x,y)*vy(x,y)* dyvy(x,y)* px(x,y) - 7*wx(x,y)*wy(x,y)*vy(x,y)*dyvy(x,y)* px(x,y) - 9*A*wx(x,y)*wy(x,y)*vy(x,y)* dyvy(x,y)* px(x,y) + 2*wy(x,y)^2*vy(x,y)*dyvy(x,y)* px(x,y) + 2*A*wy(x,y)^2*vy(x,y)* dyvy(x,y)* px(x,y) - 13*wy(x,y)*vy(x,y)*dxvx(x,y)* px(x,y) - 15*A*wy(x,y)*vy(x,y)* dxvx(x,y)* px(x,y) + 13*wx(x,y)*wy(x,y)*vy(x,y)*dxvx(x,y)* px(x,y) + 15*A*wx(x,y)*wy(x,y)*vy(x,y)* dxvx(x,y)* px(x,y) - wy(x,y)^2*vy(x,y)*dxvx(x,y)* px(x,y) + A*wy(x,y)^2*vy(x,y)* dxvx(x,y)* px(x,y) - 9*wx(x,y)*vy(x,y)*dxvy(x,y)* py(x,y) - 11*A*wx(x,y)*vy(x,y)* dxvy(x,y)* py(x,y) + 9*wx(x,y)^2*vy(x,y)*dxvy(x,y)* py(x,y) + 11*A*wx(x,y)^2*vy(x,y)* dxvy(x,y)* py(x,y) - 5*wx(x,y)*wy(x,y)*vy(x,y)*dxvy(x,y)* py(x,y) - 3*A*wx(x,y)*wy(x,y)*vy(x,y)* dxvy(x,y)* py(x,y) - 5*wx(x,y)*px(x,y)* py(x,y) - 3*A*wx(x,y)*px(x,y)*py(x,y) + 5*wx(x,y)^2*px(x,y)* py(x,y) + 3*A*wx(x,y)^2*px(x,y)* py(x,y) - 5*wy(x,y)*px(x,y)*py(x,y) - 3*A*wy(x,y)*px(x,y)* py(x,y) - 4*wx(x,y)*wy(x,y)*px(x,y)* py(x,y) - 12*A*wx(x,y)*wy(x,y)* px(x,y)* py(x,y) + 4*A*wx(x,y)^2*wy(x,y)*px(x,y)* py(x,y) + 5*wy(x,y)^2*px(x,y)*py(x,y) + 3*A*wy(x,y)^2*px(x,y)* py(x,y) + 4*A*wx(x,y)*wy(x,y)^2*px(x,y)* py(x,y) - dyvx(x,y)*(wx(x,y)*((-(5 + 3*A))*(-1 + wx(x,y)) + (3 + 5*A + 6*(1 + A)*wx(x,y))* wy(x,y) - 8*(1 + A)*wy(x,y)^2)* dyvy(x,y) + wy(x,y)*((-(-1 + wx(x,y)))*(9 + 11*A + 4*(1 + A)*wx(x,y)) + (-9 - 11*A + 10*(1 + A)*wx(x,y))*wy(x,y))* dxvx(x,y) + ((-(5 + 3*A))*(-1 + wx(x,y))* wx(x,y) + (2 + 2*A + 7*wx(x,y) + 9*A*wx(x,y))*wy(x,y) - 2*(1 + A)*wy(x,y)^2)*vy(x,y)* py(x,y)) + 2*wy(x,y)*dypx(x,y) + 2*A*wy(x,y)*dypx(x,y) - 2*wx(x,y)^2*wy(x,y)*dypx(x,y) - 2*A*wx(x,y)^2*wy(x,y)*dypx(x,y) - 4*wy(x,y)^2*dypx(x,y) - 4*A*wy(x,y)^2*dypx(x,y) + 2*wy(x,y)^3*dypx(x,y) + 2*A*wy(x,y)^3*dypx(x,y) + 2*wx(x,y)*dxpy(x,y) + 2*A*wx(x,y)*dxpy(x,y) - 4*wx(x,y)^2*dxpy(x,y) - 4*A*wx(x,y)^2*dxpy(x,y) + 2*wx(x,y)^3*dxpy(x,y) + 2*A*wx(x,y)^3*dxpy(x,y) - 2*wx(x,y)*wy(x,y)^2*dxpy(x,y) - 2*A*wx(x,y)*wy(x,y)^2*dxpy(x,y) - 2*(1 + A)*wy(x,y)*(-1 + wx(x,y) + wy(x,y))*vy(x,y)* dtdtvx(x,y) + vx(x,y)*(2*(1 + A)*wx(x,y)*(-1 + wy(x,y))*(-1 + wx(x,y) + wy(x,y))* dydyvy(x,y) + 2*wy(x,y)*dxdyvx(x,y) + 2*A*wy(x,y)*dxdyvx(x,y) - 2*wx(x,y)^2*wy(x,y)*dxdyvx(x,y) - 2*A*wx(x,y)^2*wy(x,y)* dxdyvx(x,y) - 4*wy(x,y)^2* dxdyvx(x,y) - 4*A*wy(x,y)^2* dxdyvx(x,y) + 2*wy(x,y)^3* dxdyvx(x,y) + 2*A*wy(x,y)^3* dxdyvx(x,y) + 2*wx(x,y)*wy(x,y)* dxdxvy(x,y) + 2*A*wx(x,y)*wy(x,y)* dxdxvy(x,y) - 2*wx(x,y)^2*wy(x,y)* dxdxvy(x,y) - 2*A*wx(x,y)^2*wy(x,y)* dxdxvy(x,y) - 2*wx(x,y)*wy(x,y)^2* dxdxvy(x,y) - 2*A*wx(x,y)*wy(x,y)^2* dxdxvy(x,y) - 9*wy(x,y)*dyvx(x,y)*px(x,y) - 11*A*wy(x,y)*dyvx(x,y)* px(x,y) - 5*wx(x,y)*wy(x,y)*dyvx(x,y)* px(x,y) - 3*A*wx(x,y)*wy(x,y)* dyvx(x,y)* px(x,y) + 9*wy(x,y)^2*dyvx(x,y)* px(x,y) + 11*A*wy(x,y)^2* dyvx(x,y)* px(x,y) - 2*wx(x,y)*dxvy(x,y)* px(x,y) - 2*A*wx(x,y)* dxvy(x,y)* px(x,y) + 2*wx(x,y)^2*dxvy(x,y)* px(x,y) + 2*A*wx(x,y)^2* dxvy(x,y)* px(x,y) - 5*wy(x,y)*dxvy(x,y)* px(x,y) - 3*A*wy(x,y)* dxvy(x,y)* px(x,y) - 7*wx(x,y)*wy(x,y)*dxvy(x,y)* px(x,y) - 9*A*wx(x,y)*wy(x,y)* dxvy(x,y)* px(x,y) + 5*wy(x,y)^2*dxvy(x,y)* px(x,y) + 3*A*wy(x,y)^2* dxvy(x,y)* px(x,y) - 13*wx(x,y)*dyvy(x,y)* py(x,y) - 15*A*wx(x,y)* dyvy(x,y)* py(x,y) - wx(x,y)^2*dyvy(x,y)* py(x,y) + A*wx(x,y)^2* dyvy(x,y)* py(x,y) + 13*wx(x,y)*wy(x,y)*dyvy(x,y)* py(x,y) + 15*A*wx(x,y)*wy(x,y)* dyvy(x,y)* py(x,y) - 2*wx(x,y)*dxvx(x,y)* py(x,y) - 2*A*wx(x,y)* dxvx(x,y)* py(x,y) + 2*wx(x,y)^2*dxvx(x,y)* py(x,y) + 2*A*wx(x,y)^2* dxvx(x,y)* py(x,y) - 5*wy(x,y)*dxvx(x,y)* py(x,y) - 3*A*wy(x,y)* dxvx(x,y)* py(x,y) - 7*wx(x,y)*wy(x,y)*dxvx(x,y)* py(x,y) - 9*A*wx(x,y)*wy(x,y)* dxvx(x,y)* py(x,y) + 5*wy(x,y)^2*dxvx(x,y)* py(x,y) + 3*A*wy(x,y)^2* dxvx(x,y)* py(x,y) + vy(x,y)*((-(1 + A))*(2 + 5*wx(x,y) - 2*wy(x,y))*wy(x,y)* dyvx(x,y)^2 - (1 + A)* wx(x,y)*(6 + wx(x,y) - 6*wy(x,y))* dyvy(x,y)^2 + 2*(1 + A)*( (-1 + wx(x,y))*wx(x,y) - (1 + 5*wx(x,y))*wy(x,y) + wy(x,y)^2)* dyvy(x,y)* dxvx(x,y) - 6*wy(x,y)*dxvx(x,y)^2 - 6*A*wy(x,y)*dxvx(x,y)^2 + 6*wx(x,y)*wy(x,y)*dxvx(x,y)^2 + 6*A*wx(x,y)*wy(x,y)*dxvx(x,y)^2 - wy(x,y)^2*dxvx(x,y)^2 - A*wy(x,y)^2*dxvx(x,y)^2 + 2*(1 + A)*(wx(x,y)^2 + (-1 + wy(x,y))*wy(x,y) - wx(x,y)*(1 + 5*wy(x,y)))* dyvx(x,y)* dxvy(x,y) - 2*wx(x,y)*dxvy(x,y)^2 - 2*A*wx(x,y)*dxvy(x,y)^2 + 2*wx(x,y)^2*dxvy(x,y)^2 + 2*A*wx(x,y)^2*dxvy(x,y)^2 - 5*wx(x,y)*wy(x,y)*dxvy(x,y)^2 - 5*A*wx(x,y)*wy(x,y)*dxvy(x,y)^2 - 7*wy(x,y)*px(x,y)^2 - 9*A*wy(x,y)*px(x,y)^2 + 2*A*wx(x,y)*wy(x,y)* px(x,y)^2 + 2*A*wy(x,y)^2* px(x,y)^2 - 7*wx(x,y)*py(x,y)^2 - 9*A*wx(x,y)*py(x,y)^ 2 + 2*A*wx(x,y)^2*py(x,y)^2 + 2*A*wx(x,y)*wy(x,y)*py(x,y)^2 + 4*wx(x,y)*dypy(x,y) + 4*A*wx(x,y)*dypy(x,y) - 4*wx(x,y)^2*dypy(x,y) - 4*A*wx(x,y)^2* dypy(x,y) - 4*wx(x,y)*wy(x,y)*dypy(x,y) - 4*A*wx(x,y)*wy(x,y)*dypy(x,y) - 4*(1 + A)*wy(x,y)*(-1 + wx(x,y) + wy(x,y))* dxpx(x,y)) - 2*(1 + A)*wx(x,y)*(-1 + wx(x,y) + wy(x,y))* dtdtvy(x,y)))));

      B(x,y) = -((1/(9*wx(x,y)*wy(x,y)*(-1 + wx(x,y) + wy(x,y))^(7/2)))* (pi* T(x,y)*(-6*wy(x,y)* sqrt(-1 + wx(x,y) + wy(x,y))*(1 + A + A*wx(x,y) - (1 + A)*wy(x,y))* dyvx(x,y)^2 - 8*pi*(-1 + wx(x,y) + wy(x,y))^2* T(x,y)*(wx(x,y)*dyvy(x,y) - wy(x,y)*(dxvx(x,y) + vx(x,y)* px(x,y)) + wx(x,y)*vy(x,y)* py(x,y)) - 3*sqrt(-1 + wx(x,y) + wy(x,y))* dyvx(x,y)*(2*(1 + A)*(wx(x,y) - wy(x,y))* (-1 + wx(x,y) + wy(x,y))* dxvy(x,y) + (9 + 11*A - 5*wx(x,y) - 3*A*wx(x,y))*wy(x,y)*vy(x,y)* px(x,y) + vx(x,y)*(vy(x,y)*((-(5 + 3*A))*wx(x,y)* dyvy(x,y) + (9 + 11*A)*wy(x,y)* dxvx(x,y)) + (-2*(1 + A)*(-1 + wy(x,y))* wy(x,y) + wx(x,y)*(-5 - 3*A + (3 + A)*wy(x,y)))* py(x,y))) + sqrt(-1 + wx(x,y) + wy(x,y))*(3*(1 + A)*wx(x,y)*(6 + wx(x,y) - 6*wy(x,y))* dyvy(x,y)^2 - 6*(1 + A)*wy(x,y)^2*vy(x,y)* dxdyvx(x,y) + 6*(1 + A)*wx(x,y)^2*vx(x,y)* dxdyvy(x,y) - 3*dyvy(x,y)*(2*(1 + A)*(wx(x,y) - wy(x,y))*(-1 + wx(x,y) + wy(x,y))*dxvx(x,y) + vx(x,y)*((-(9 + 11*A))*wx(x,y)*vy(x,y)* dxvy(x,y) + (-2*(1 + A)*(-1 + wy(x,y))*wy(x,y) + wx(x,y)*(-5 - 3*A + (3 + A)*wy(x,y)))* px(x,y)) + wx(x,y)*(-13 - wx(x,y) + 4*wy(x,y) + A*(-15 + wx(x,y) + 4*wy(x,y)))*vy(x,y)* py(x,y)) + wx(x,y)*vx(x,y)* (-6*(1 + A)*dxdyvy(x,y) + 3*(9 + 11*A)* dxvy(x,y)* py(x,y)) + 3*wy(x,y)*vy(x,y)*(2*(1 + A)*dxdyvx(x,y) - (5 + 3*A)*(dxvy(x,y)* px(x,y) + dxvx(x,y)* py(x,y))) - 6*(1 + A)*wy(x,y)^2*vx(x,y)*vy(x,y)* dypx(x,y) - 3*wy(x,y)*vx(x,y)*vy(x,y)* ((5 + 3*A)*(dxvx(x,y)* dxvy(x,y) + px(x,y)* py(x,y)) - 2*(1 + A)* dypx(x,y)) + 3*wx(x,y)*(2*(1 + A)*dxvy(x,y)^2 + (7 + 9*A)* py(x,y)^2 - 4*(1 + A)*dypy(x,y)) - 6*wx(x,y)^2*((1 + A)*dxvy(x,y)^2 + A*py(x,y)^2 - 2*(1 + A)* dypy(x,y)) + 3*wx(x,y)*wy(x,y)*(6*(1 + A)*dxvx(x,y)^2 + 2*A*dxvy(x,y)^2 + (1 + A)*(5*px(x,y)^2 - 5* py(x,y)^2 + 6*dypy(x,y) - 6*dxpx(x,y))) - 6*(1 + A)*wx(x,y)^2*wy(x,y)*(dypy(x,y) - dxpx(x,y)) - 6*(1 + A)*wx(x,y)*wy(x,y)^2* (dypy(x,y) - dxpx(x,y)) - 3*wy(x,y)*(6*(1 + A)*dxvx(x,y)^2 + (7 + 9*A)* px(x,y)^2 - 4*(1 + A)*dxpx(x,y)) - 3*wy(x,y)^2*((1 + A)*dxvx(x,y)^2 - 2*A*px(x,y)^2 + 4*(1 + A)* dxpx(x,y)) - 6*(1 + A)*wx(x,y)*wy(x,y)*vx(x,y)* vy(x,y)*(dypx(x,y) - dxpy(x,y)) + 6*(1 + A)*wx(x,y)^2*vx(x,y)* vy(x,y)* dxpy(x,y) + wx(x,y)*vx(x,y)* vy(x,y)*(3*(5 + 3*A)*px(x,y)* py(x,y) - 6*(1 + A)* dxpy(x,y)) + 3*wx(x,y)*wy(x,y)*vx(x,y)* (2*(1 + A)*dxdyvy(x,y) - 2*(1 + A)*dxdxvx(x,y) + 4*dxvx(x,y)* px(x,y) + 4*A*dxvx(x,y)* px(x,y) - 5*dxvy(x,y)* py(x,y) - 3*A*dxvy(x,y)* py(x,y) - 2*(1 + A)*dtdtvx(x,y)) - 3*wy(x,y)^2*vx(x,y)*(2*(1 + A)*dxdxvx(x,y) - (-1 + A)*dxvx(x,y)* px(x,y) + 2*(1 + A)*dtdtvx(x,y)) + wy(x,y)*vx(x,y)*(6*(1 + A)*dxdxvx(x,y) - 3*(13 + 15*A)*dxvx(x,y)* px(x,y) + 6*(1 + A)* dtdtvx(x,y)) - 6*(1 + A)*wx(x,y)*vy(x,y)* (dydyvy(x,y) - dxvy(x,y)*px(x,y) - dxvx(x,y)* py(x,y) + dtdtvy(x,y)) + 6*(1 + A)*wx(x,y)^2*vy(x,y)* (dydyvy(x,y) - dxvy(x,y)*px(x,y) - dxvx(x,y)* py(x,y) + dtdtvy(x,y)) + 3*wx(x,y)*wy(x,y)*vy(x,y)* (2*(1 + A)*dydyvy(x,y) - 2*(1 + A)*dxdyvx(x,y) + (3 + A)*(dxvy(x,y)* px(x,y) + dxvx(x,y)* py(x,y)) + 2*(1 + A)* dtdtvy(x,y))))));

    end
  end


  dxPi = dx(Pi);
  dyPi = dy(Pi);
  dxB = dx(B);
  dyB = dy(B);

  end

  if (ZERO_HYDRO)
    T = T + DT * dtT;
    vx = vx + DT * px;
    vy = vy + DT * py;
    phydro()
    pause(0.003);
    sumT(t)=sum(sum(T));
  end
  end

  toc

  if (ZERO_HYDRO)
    r = sumT;
    return;
  end

  for t = 1:TIME_STEPS
    disp('time step ');
    disp(t);
    for x = 1:GRID_SIZE
      for y = 1:GRID_SIZE

        dtpx(x,y) = (-3*B(x,y)* vx(x,y)*(-1 + vx(x,y)^2)*(-1 + vy(x,y)^2)*(-1 + vx(x,y)^2 + vy(x,y)^2)^3 + 3*vy(x,y)^9*Pi(x,y) + 2*(1 + A)*pi*T(x,y)* vy(x,y)^8*(vx(x,y)*dxvx(x,y) + px(x,y)) + 2*vy(x,y)^7*(6*(-1 + vx(x,y)^2)*Pi(x,y) - (1 + A)*pi* T(x,y)*((-1 + vx(x,y)^2)* dyvx(x,y) + (-1 + vx(x,y)^2)*dxvy(x,y) + vx(x,y)*dxdyvx(x,y) + dypx(x,y))) + pi*T(x,y)* vy(x,y)^6*(-2*A*dydyvx(x,y) + 2*(1 + A)* vx(x,y)^3*(dyvy(x,y) + 2*dxvx(x,y)) - 2*A*dxdyvy(x,y) - 4*A*px(x,y) + 2*dyvy(x,y)* px(x,y) + 2*A*dyvy(x,y)* px(x,y) - dxvx(x,y)* px(x,y) + A*dxvx(x,y)* px(x,y) - 2*(dydyvx(x,y) + dxdyvy(x,y) + 2*px(x,y)) + 2*(1 + A)* vx(x,y)^2*(dydyvx(x,y) + dxdyvy(x,y) - dxdxvx(x,y) + 3*px(x,y)) + 2*(1 + A)*dyvx(x,y)* py(x,y) + vx(x,y)*(2*(1 + A)*dyvx(x,y)^2 + 2*(1 + A)* dyvy(x,y)*(-1 + dxvx(x,y)) - 2*dxvx(x,y) - 2*A*dxvx(x,y) - dxvx(x,y)^2 - A*dxvx(x,y)^2 + 2*(1 + A)*dyvx(x,y)* dxvy(x,y) + 2*A*px(x,y)^2 - 4*(1 + A)*dxpx(x,y))) - pi*T(x,y)* vx(x,y)*(-2*A*(-1 + vx(x,y)^2)* dyvx(x,y)^2 + 2*(1 + A)*vx(x,y)^6*dyvy(x,y) - dyvy(x,y)^2 - A*dyvy(x,y)^2 + dxvx(x,y)^2 + A*dxvx(x,y)^2 - 2*A*dxvy(x,y)^2 - 2*(1 + A)* vx(x,y)^5*(dxdyvy(x,y) - px(x,y)) + 2*px(x,y)^2 + 2*A*px(x,y)^2 - 2*py(x,y)^2 - 2*A*py(x,y)^2 + 2*(1 + A)*vx(x,y)*(-1 + vx(x,y)^2)* dyvx(x,y)*(vx(x,y)*dxvy(x,y) + py(x,y)) + 2* vx(x,y)*((1 + A)*(-dxdyvy(x,y) + dxdxvx(x,y) - (dyvy(x,y) - 5*dxvx(x,y))* px(x,y)) - 2*(1 + 2*A)*dxvy(x,y)* py(x,y)) + 2* vx(x,y)^3*((1 + A)*(2*dxdyvy(x,y) - dxdxvx(x,y) + (-1 + dyvy(x,y) + 2*dxvx(x,y))* px(x,y)) + 2*(1 + 2*A)*dxvy(x,y)* py(x,y)) - 2*dypy(x,y) - 2*A*dypy(x,y) + vx(x,y)^2*(2*(1 + A)*dyvy(x,y)^2 - 2*(1 + A)*dyvy(x,y)* (-1 + dxvx(x,y)) - 2*(1 + A)*dxvx(x,y) + 6*(1 + A)*dxvx(x,y)^2 - 2*dxvy(x,y)^2 + 5*(1 + A)*px(x,y)^2 + 4*A*dypy(x,y) + 2*(py(x,y)^2 + 2*dypy(x,y))) + 2*dxpx(x,y) + 2*A*dxpx(x,y) - vx(x,y)^4*((1 + A)*dyvy(x,y)^2 - 2*(1 + A)* dyvy(x,y)*(-2 + dxvx(x,y)) - 2*(1 + A)*(dxvx(x,y) + dxvy(x,y)^2) - 2*A*py(x,y)^2 + 2*(1 + A)*dypy(x,y) + 2*(1 + A)*dxpx(x,y))) + pi*T(x,y)* vy(x,y)^4*(4*A*dydyvx(x,y) + 2*(1 + A)* vx(x,y)^5*(2*dyvy(x,y) + dxvx(x,y)) + 4*A*dxdyvy(x,y) + 4*(dydyvx(x,y) + dxdyvy(x,y)) + 2*px(x,y) + 2*A*px(x,y) + 5*dyvy(x,y)* px(x,y) + 7*A*dyvy(x,y)* px(x,y) + 2*dxvx(x,y)* px(x,y) - 2*A*dxvx(x,y)* px(x,y) + 2*(1 + A)* vx(x,y)^4*(dydyvx(x,y) - dxdxvx(x,y) + 3*px(x,y)) + ((5 + 7*A)* dyvx(x,y) + (5 + 3*A)* dxvy(x,y))* py(x,y) + vx(x,y)^2*(-6*(1 + A)*dydyvx(x,y) - 4*(1 + A)*dxdyvy(x,y) + 2*dxdxvx(x,y) + 2*A*dxdxvx(x,y) - 10*px(x,y) - 10*A*px(x,y) - 5*dyvy(x,y)* px(x,y) - 7*A*dyvy(x,y)* px(x,y) + 12*dxvx(x,y)* px(x,y) + 16*A*dxvx(x,y)* px(x,y) - ((5 + 7*A)*dyvx(x,y) + (5 + 3*A)* dxvy(x,y))* py(x,y)) - vx(x,y)^3*(5*(1 + A)*dyvx(x,y)^2 - 6*(1 + A)*dyvy(x,y)^2 + 2*dxvx(x,y) + 2*A*dxvx(x,y) - 6*dxvx(x,y)^2 - 6*A*dxvx(x,y)^2 + 2*(1 + A)* dyvy(x,y)*(4 + 5*dxvx(x,y)) + 10*(1 + A)*dyvx(x,y)* dxvy(x,y) + 5*dxvy(x,y)^2 + 5*A*dxvy(x,y)^2 - 2*A*px(x,y)^2 - 2*A*py(x,y)^2 + 2*dypy(x,y) + 2*A*dypy(x,y) + 6*(1 + A)*dxpx(x,y)) + vx(x,y)*(5*(1 + A)*dyvx(x,y)^2 - 6*(1 + A)*dyvy(x,y)^2 - 2*dxvx(x,y) - 2*A*dxvx(x,y) + dxvx(x,y)^2 + A*dxvx(x,y)^2 + 2*(1 + A)*dyvy(x,y)* (2 + 5*dxvx(x,y)) + 10*(1 + A)*dyvx(x,y)* dxvy(x,y) + 5*dxvy(x,y)^2 + 5*A*dxvy(x,y)^2 + 7*px(x,y)^2 + 5*A*px(x,y)^2 - 2*A*py(x,y)^2 + 2*dypy(x,y) + 2*A*dypy(x,y) + 6*(1 + A)*dxpx(x,y))) + pi*T(x,y)* vy(x,y)^2*(2*(1 + A)*vx(x,y)^7* dyvy(x,y) - 2*A*dydyvx(x,y) - 2*A*dxdyvy(x,y) - 2*(dydyvx(x,y) + dxdyvy(x,y)) - 2*(1 + A)* vx(x,y)^6*(dxdyvy(x,y) - px(x,y)) - 7*dyvy(x,y)* px(x,y) - 9*A*dyvy(x,y)* px(x,y) - dxvx(x,y)* px(x,y) + A*dxvx(x,y)* px(x,y) - ((7 + 9*A)*dyvx(x,y) + (5 + 3*A)* dxvy(x,y))* py(x,y) + vx(x,y)^2*(4*(1 + A)*dydyvx(x,y) + 2*(1 + A)*dxdxvx(x,y) + (4*(1 + A) + (3 + 5*A)* dyvy(x,y) - 2*(1 + 3*A)*dxvx(x,y))* px(x,y) + ((3 + 5*A)*dyvx(x,y) + (1 - 5*A)* dxvy(x,y))* py(x,y)) + 2* vx(x,y)^4*((1 + A)*(-dydyvx(x,y) + 2*dxdyvy(x,y) + (-4 + dyvy(x,y) + 2*dxvx(x,y))* px(x,y)) + ((1 + A)* dyvx(x,y) + 2*(1 + 2*A)*dxvy(x,y))* py(x,y)) + vx(x,y)*((-(7 + 5*A))*dyvx(x,y)^2 + 5*(1 + A)*dyvy(x,y)^2 + 2*dxvx(x,y) + 2*A*dxvx(x,y) + dxvx(x,y)^2 + A*dxvx(x,y)^2 - 2*(1 + A)*dyvy(x,y)* (1 + 6*dxvx(x,y)) - 12*(1 + A)*dyvx(x,y)* dxvy(x,y) - 5*dxvy(x,y)^2 - 7*A*dxvy(x,y)^2 - 5*px(x,y)^2 - 5*A*px(x,y)^2 - 2*py(x,y)^2 - 4*(1 + A)*dypy(x,y)) - vx(x,y)^5*((1 + A)*dyvy(x,y)^2 - 2*(1 + A)*dyvy(x,y)* (-4 + dxvx(x,y)) - 2*(1 + A)* dxvy(x,y)*(dyvx(x,y) + dxvy(x,y)) - 2*A*py(x,y)^2 + 2*(1 + A)*dypy(x,y) + 2*(1 + A)*dxpx(x,y)) + vx(x,y)^3*((5 + 3*A)*dyvx(x,y)^2 - 4*(1 + A)*dyvy(x,y)^2 - 4*dxvx(x,y) - 4*A*dxvx(x,y) + 8*(1 + A)* dyvy(x,y)*(1 + dxvx(x,y)) + 8*(1 + A)*dyvx(x,y)* dxvy(x,y) + 3*dxvy(x,y)^2 + 5*A*dxvy(x,y)^2 + 5*px(x,y)^2 + 3*A*px(x,y)^2 + 2*py(x,y)^2 - 2*A*py(x,y)^2 + 6*dypy(x,y) + 6*A*dypy(x,y) + 6*(1 + A)*dxpx(x,y))) + vy(x,y)^5*(18*(-1 + vx(x,y)^2)^2*Pi(x,y) + pi*T(x,y)*(-4*dxvy(x,y) - 4*A*dxvy(x,y) - 4*(1 + A)*vx(x,y)^4*dxvy(x,y) + 4*dyvy(x,y)* dxvy(x,y) + 4*A*dyvy(x,y)* dxvy(x,y) - dxvx(x,y)* dxvy(x,y) - 3*A*dxvx(x,y)* dxvy(x,y) - 2*(1 + A)* vx(x,y)^3*(dydyvy(x,y) - dxdxvy(x,y)) + dyvx(x,y)*(-4*(1 + A)*vx(x,y)^4 + 4*(1 + A)*(-1 + 2*dyvy(x,y)) - 2*(1 + A)* vx(x,y)^2*(-4 + 4*dyvy(x,y) - 5*dxvx(x,y)) + (-1 + A)* dxvx(x,y) + (9 + 11*A)*vx(x,y)* px(x,y)) + 5*px(x,y)* py(x,y) + 7*A*px(x,y)* py(x,y) + vx(x,y)*(2*(1 + A)*dydyvy(x,y) - 2*(1 + A)*dxdxvy(x,y) + (5 + 3*A)*(dxvy(x,y)* px(x,y) + dxvx(x,y)* py(x,y))) + 2*dypx(x,y) + 2*A*dypx(x,y) - 2*vx(x,y)^2*((1 + A)*(-4 + 2*dyvy(x,y) - 3*dxvx(x,y))* dxvy(x,y) + 2*A*px(x,y)* py(x,y) + (1 + A)*(dypx(x,y) - dxpy(x,y))) - 2*dxpy(x,y) - 2*A*dxpy(x,y))) + vy(x,y)^3*(12*(-1 + vx(x,y)^2)^3*Pi(x,y) + pi*T(x,y)*((-1 + vx(x,y)^2)*dyvx(x,y)* (-2*(1 + A) - 2*(1 + A)*vx(x,y)^4 + (7 + 5*A)* dyvy(x,y) + 2*(1 + A)* vx(x,y)^2*(3 + 3*dyvy(x,y) - 2*dxvx(x,y)) + 2*(-1 + A)*dxvx(x,y)) + 2*dxvy(x,y) + 2*A*dxvy(x,y) - 2*(1 + A)*vx(x,y)^6*dxvy(x,y) - 3*dyvy(x,y)* dxvy(x,y) - 5*A*dyvy(x,y)* dxvy(x,y) + 2*dxvx(x,y)* dxvy(x,y) + 6*A*dxvx(x,y)* dxvy(x,y) - 2*(1 + A)* vx(x,y)^5*(dydyvy(x,y) - dxdyvx(x,y) - dxdxvy(x,y)) - px(x,y)* py(x,y) - 3*A*px(x,y)* py(x,y) + vx(x,y)^3*(6*(1 + A)*dydyvy(x,y) - 6*(1 + A)*(dxdyvx(x,y) + dxdxvy(x,y)) - 4*(1 + 2*A)*dxvy(x,y)* px(x,y) + ((9 + 11*A)*dyvy(x,y) - 4*(1 + 2*A)*dxvx(x,y))* py(x,y)) + vx(x,y)*(-2*(1 + A)*(2*dydyvy(x,y) - 3*dxdyvx(x,y) - 2*dxdxvy(x,y)) + 4*(1 + 2*A)*dxvy(x,y)* px(x,y) + ((-(9 + 11*A))*dyvy(x,y) + 4*(1 + 2*A)*dxvx(x,y))* py(x,y)) + 2*dypx(x,y) + 2*A*dypx(x,y) + 4*(1 + A)*dxpy(x,y) + vx(x,y)^2*((-(8*(1 + A) + (7 + 5*A)* dyvy(x,y) - 2*(3 + A)*dxvx(x,y)))* dxvy(x,y) + (1 + 7*A)*px(x,y)* py(x,y) - 6*(1 + A)*dxpy(x,y)) + 2*vx(x,y)^4*((1 + A)*(4 + 5*dyvy(x,y) - 4*dxvx(x,y))* dxvy(x,y) - 2*A*px(x,y)* py(x,y) + (1 + A)* dxpy(x,y)))) + vy(x,y)*(3*(-1 + vx(x,y)^2)^4*Pi(x,y) + pi*T(x,y)*(dyvx(x,y)*((-(1 + 3*A))*dyvy(x,y) + (-1 + A)*dxvx(x,y)) - dyvy(x,y)* dxvy(x,y) + A*dyvy(x,y)* dxvy(x,y) - dxvx(x,y)* dxvy(x,y) - 3*A*dxvx(x,y)* dxvy(x,y) + 2*(1 + A)* vx(x,y)^6*(dyvx(x,y) + dxvy(x,y)) + 2*(1 + A)* vx(x,y)^5*(dydyvy(x,y) - dxdyvx(x,y) - dxdxvy(x,y)) - 4*px(x,y)* py(x,y) - 4*A*px(x,y)* py(x,y) + vx(x,y)*(2*(1 + A)*dydyvy(x,y) - 4*(1 + A)*dxdyvx(x,y) - 2*dxdxvy(x,y) - 2*A*dxdxvy(x,y) - 9*dyvx(x,y)* px(x,y) - 11*A*dyvx(x,y)* px(x,y) - 9*dxvy(x,y)* px(x,y) - 11*A*dxvy(x,y)* px(x,y) + (9 + 11*A)*(dyvy(x,y) - dxvx(x,y))* py(x,y)) + vx(x,y)^3*(-2*(1 + A)*(2*dydyvy(x,y) - 3*dxdyvx(x,y) - 2*dxdxvy(x,y)) + 4*(1 + 2*A)*dxvy(x,y)* px(x,y) + ((-(9 + 11*A))*dyvy(x,y) + 4*(1 + 2*A)*dxvx(x,y))* py(x,y)) - 2*A*dypx(x,y) - 2*A*dxpy(x,y) - 2*(dypx(x,y) + dxpy(x,y)) - 2*vx(x,y)^4*((1 + A)* dyvx(x,y)*(2 + 3*dyvy(x,y) - 2*dxvx(x,y)) + (1 + A)*(2 + 5*dyvy(x,y) - 4*dxvx(x,y))* dxvy(x,y) - 2*A*px(x,y)* py(x,y) + (1 + A)* dxpy(x,y)) + vx(x,y)^2*(dyvx(x,y)*((7 + 9*A)*dyvy(x,y) + 2*(1 + A - 2*(3 + 4*A)*dxvx(x,y))) + ((11 + 9*A)*dyvy(x,y) + 2*(1 + A - 2*(3 + 2*A)*dxvx(x,y)))* dxvy(x,y) - px(x,y)* py(x,y) - 3*A*px(x,y)* py(x,y) + 2*dypx(x,y) + 2*A*dypx(x,y) + 4*(1 + A)*dxpy(x,y))))) * (2*sqrt(1 - vx(x,y)^2 - vy(x,y)^2))/(3*(1 + A)^2*(-1 + vy(x,y))*(1 + vy(x,y))*(-1 + vx(x,y)^2 + vy(x,y)^2)*(vx(x,y)^2 + vy(x,y)^2)^2);







        dtpy(x, y) = ((9*wx(x,y)*wy(x,y)*(-1 + wx(x,y) + wy(x,y))^(7/2)*B(x,y)*vy(x,y) + pi*wx(x,y)*T(x,y)* (6*(1 + A)*G(x,y)*(-2 + wx(x,y))*(-1 + wy(x,y))*(-1 + wx(x,y) + wy(x,y))* dydyvy(x,y) + 6*G(x,y)*wx(x,y)*wy(x,y)*dxdyvx(x,y) + 6*A*G(x,y)*wx(x,y)*wy(x,y)*dxdyvx(x,y) - 6*G(x,y)*wx(x,y)^2*wy(x,y)*dxdyvx(x,y) - 6*A*G(x,y)*wx(x,y)^2*wy(x,y)*dxdyvx(x,y) - 6*G(x,y)*wy(x,y)^2*dxdyvx(x,y) - 6*A*G(x,y)*wy(x,y)^2*dxdyvx(x,y) + 6*G(x,y)*wy(x,y)^3*dxdyvx(x,y) + 6*A*G(x,y)*wy(x,y)^3*dxdyvx(x,y) - 6*G(x,y)*wy(x,y)*dxdxvy(x,y) - 6*A*G(x,y)*wy(x,y)*dxdxvy(x,y) + 12*G(x,y)*wx(x,y)*wy(x,y)*dxdxvy(x,y) + 12*A*G(x,y)*wx(x,y)*wy(x,y)*dxdxvy(x,y) - 6*G(x,y)*wx(x,y)^2*wy(x,y)*dxdxvy(x,y) - 6*A*G(x,y)*wx(x,y)^2*wy(x,y)*dxdxvy(x,y) + 6*G(x,y)*wy(x,y)^2*dxdxvy(x,y) + 6*A*G(x,y)*wy(x,y)^2*dxdxvy(x,y) - 6*G(x,y)*wx(x,y)*wy(x,y)^2*dxdxvy(x,y) - 6*A*G(x,y)*wx(x,y)*wy(x,y)^2*dxdxvy(x,y) + 3*G(x,y)*wy(x,y)*dyvx(x,y)* px(x,y) - 15*A*G(x,y)*wy(x,y)*dyvx(x,y)* px(x,y) - 15*G(x,y)*wx(x,y)*wy(x,y)*dyvx(x,y)* px(x,y) - 9*A*G(x,y)*wx(x,y)*wy(x,y)*dyvx(x,y)* px(x,y) + 12*G(x,y)*wy(x,y)^2*dyvx(x,y)* px(x,y) + 24*A*G(x,y)*wy(x,y)^2*dyvx(x,y)* px(x,y) + 12*G(x,y)*dxvy(x,y)* px(x,y) + 12*A*G(x,y)*dxvy(x,y)* px(x,y) - 18*G(x,y)*wx(x,y)*dxvy(x,y)* px(x,y) - 18*A*G(x,y)*wx(x,y)*dxvy(x,y)* px(x,y) + 6*G(x,y)*wx(x,y)^2*dxvy(x,y)* px(x,y) + 6*A*G(x,y)*wx(x,y)^2*dxvy(x,y)* px(x,y) + 9*G(x,y)*wy(x,y)*dxvy(x,y)* px(x,y) + 15*A*G(x,y)*wy(x,y)*dxvy(x,y)* px(x,y) - 15*G(x,y)*wx(x,y)*wy(x,y)*dxvy(x,y)* px(x,y) - 21*A*G(x,y)*wx(x,y)*wy(x,y)*dxvy(x,y)* px(x,y) + 6*G(x,y)*wy(x,y)^2*dxvy(x,y)* px(x,y) + 6*A*G(x,y)*wy(x,y)^2*dxvy(x,y)* px(x,y) + (8*pi*(-2 + wx(x,y) + wy(x,y))*(-1 + wx(x,y) + wy(x,y))^2*T(x,y) + 3*G(x,y)*(((-13 + A*(-15 + wx(x,y)) - wx(x,y))*(-2 + wx(x,y)) + 2* (-15 + 6*wx(x,y) + A*(-17 + 8*wx(x,y)))*wy(x,y) + 4*(1 + A)*wy(x,y)^2)*dyvy(x,y) + (2*(1 + A)*(-2 + wx(x,y))*(-1 + wx(x,y)) + (3 + A*(5 - 7*wx(x,y)) - 5*wx(x,y))*wy(x,y) + 2*(1 + A)*wy(x,y)^2)* dxvx(x,y)))* py(x,y) + vy(x,y)*(3*G(x,y)*wy(x,y)*(3 + A - 5*wx(x,y) - 5*A*wx(x,y) + 2*(1 + A)*wy(x,y))* dyvx(x,y)^2 + 8*pi*(-1 + wx(x,y) + wy(x,y))^2* T(x,y)*((-2 + wx(x,y))*dyvy(x,y) + wy(x,y)*dxvx(x,y)) + 6*(1 + A)*G(x,y)*(2 + wx(x,y)^2 + wy(x,y)*(4 + wy(x,y)) - wx(x,y)*(3 + 5*wy(x,y)))* dyvx(x,y)* dxvy(x,y) - 3*G(x,y)*((1 + A)*(-2 + wx(x,y))*(6 + wx(x,y) - 6*wy(x,y))* dyvy(x,y)^2 - 2*(1 + A)*(2 + wx(x,y)^2 + wy(x,y)*(4 + wy(x,y)) - wx(x,y)*(3 + 5*wy(x,y)))* dyvy(x,y)* dxvx(x,y) - (1 + A)*(-6 + 6*wx(x,y) - wy(x,y))*wy(x,y)* dxvx(x,y)^2 - 4*dxvy(x,y)^2 - 4*A*dxvy(x,y)^2 + 6*wx(x,y)*dxvy(x,y)^2 + 6*A*wx(x,y)*dxvy(x,y)^2 - 2*wx(x,y)^2*dxvy(x,y)^2 - 2*A*wx(x,y)^2*dxvy(x,y)^2 - 5*wy(x,y)*dxvy(x,y)^2 - 7*A*wy(x,y)*dxvy(x,y)^2 + 5*wx(x,y)*wy(x,y)*dxvy(x,y)^2 + 5*A*wx(x,y)*wy(x,y)*dxvy(x,y)^2 + 2*wy(x,y)*px(x,y)^2 + 6*A*wy(x,y)*px(x,y)^2 - 2*A*wx(x,y)*wy(x,y)*px(x,y)^2 - 2*A*wy(x,y)^2*px(x,y)^2 - 14*py(x,y)^2 - 18*A*py(x,y)^2 + 7*wx(x,y)*py(x,y)^2 + 13*A*wx(x,y)*py(x,y)^2 - 2*A*wx(x,y)^2*py(x,y)^2 + 5*wy(x,y)*py(x,y)^2 + 7*A*wy(x,y)*py(x,y)^2 - 2*A*wx(x,y)*wy(x,y)*py(x,y)^2 + 8*dypy(x,y) + 8*A*dypy(x,y) - 12*wx(x,y)*dypy(x,y) - 12*A*wx(x,y)*dypy(x,y) + 4*wx(x,y)^2*dypy(x,y) + 4*A*wx(x,y)^2*dypy(x,y) - 10*wy(x,y)*dypy(x,y) - 10*A*wy(x,y)*dypy(x,y) + 6*wx(x,y)*wy(x,y)*dypy(x,y) + 6*A*wx(x,y)*wy(x,y)*dypy(x,y) + 2*wy(x,y)^2*dypy(x,y) + 2*A*wy(x,y)^2*dypy(x,y) + 2*(1 + A)*wy(x,y)*(-1 + wx(x,y) + wy(x,y))* dxpx(x,y)))) - vx(x,y)*(9*(-1 + wx(x,y) + wy(x,y))^(9/2)*Pi(x,y) + pi*wx(x,y)*T(x,y)* (8*pi*wy(x,y)*(-1 + wx(x,y) + wy(x,y))^2* T(x,y)*(dyvx(x,y) + dxvy(x,y)) + 3*G(x,y)*(-18*dyvy(x,y)* dxvy(x,y) - 22*A*dyvy(x,y)* dxvy(x,y) + 9*wx(x,y)*dyvy(x,y)* dxvy(x,y) + 11*A*wx(x,y)*dyvy(x,y)* dxvy(x,y) + 14*wy(x,y)*dyvy(x,y)* dxvy(x,y) + 18*A*wy(x,y)*dyvy(x,y)* dxvy(x,y) - 10*wx(x,y)*wy(x,y)*dyvy(x,y)* dxvy(x,y) - 10*A*wx(x,y)*wy(x,y)*dyvy(x,y)* dxvy(x,y) + 4*wy(x,y)^2*dyvy(x,y)* dxvy(x,y) + 4*A*wy(x,y)^2*dyvy(x,y)* dxvy(x,y) - 3*wy(x,y)*dxvx(x,y)* dxvy(x,y) - 5*A*wy(x,y)*dxvx(x,y)* dxvy(x,y) + 8*wx(x,y)*wy(x,y)*dxvx(x,y)* dxvy(x,y) + 8*A*wx(x,y)*wy(x,y)*dxvx(x,y)* dxvy(x,y) - 6*wy(x,y)^2*dxvx(x,y)* dxvy(x,y) - 6*A*wy(x,y)^2*dxvx(x,y)* dxvy(x,y) - 10*px(x,y)* py(x,y) - 6*A*px(x,y)* py(x,y) + 5*wx(x,y)*px(x,y)* py(x,y) + 3*A*wx(x,y)*px(x,y)* py(x,y) + wy(x,y)*px(x,y)* py(x,y) - 9*A*wy(x,y)*px(x,y)* py(x,y) + 4*A*wx(x,y)*wy(x,y)*px(x,y)* py(x,y) + 4*A*wy(x,y)^2*px(x,y)* py(x,y) + dyvx(x,y)*(((5 + 3*A)*(-2 + wx(x,y)) - 2*(-1 + A + 3*(1 + A)*wx(x,y))*wy(x,y) + 8*(1 + A)*wy(x,y)^2)* dyvy(x,y) + wy(x,y)*(5 + 7*A + 4*wx(x,y) + 4*A*wx(x,y) - 10*(1 + A)*wy(x,y))* dxvx(x,y) + ((5 + 3*A)*(-2 + wx(x,y)) - 4*(1 + 2*A)*wy(x,y))* vy(x,y)* py(x,y)) - vy(x,y)*(2*(1 + A)*wy(x,y)*(-1 + wx(x,y) + wy(x,y))* dydyvx(x,y) - 2*(1 + A)*(-2 + wx(x,y) - wy(x,y))*(-1 + wx(x,y) + wy(x,y))* dxdyvy(x,y) + 2*wy(x,y)*dxdxvx(x,y) + 2*A*wy(x,y)*dxdxvx(x,y) - 2*wx(x,y)* wy(x,y)*dxdxvx(x,y) - 2*A*wx(x,y)*wy(x,y)*dxdxvx(x,y) - 2*wy(x,y)^2* dxdxvx(x,y) - 2*A*wy(x,y)^2*dxdxvx(x,y) + 10* dyvy(x,y)* px(x,y) + 6*A*dyvy(x,y)* px(x,y) - 5*wx(x,y)*dyvy(x,y)* px(x,y) - 3*A*wx(x,y)*dyvy(x,y)* px(x,y) + 4*wy(x,y)*dyvy(x,y)* px(x,y) + 8*A*wy(x,y)* dyvy(x,y)* px(x,y) - 9*wy(x,y)*dxvx(x,y)* px(x,y) - 11*A*wy(x,y)*dxvx(x,y)* px(x,y) - (9 + 11*A)*(-2 + wx(x,y))* dxvy(x,y)* py(x,y)) + 2*wy(x,y)*dypx(x,y) + 2*A*wy(x,y)*dypx(x,y) - 2*wx(x,y)*wy(x,y)*dypx(x,y) - 2*A*wx(x,y)*wy(x,y)*dypx(x,y) - 2*wy(x,y)^2*dypx(x,y) - 2*A*wy(x,y)^2*dypx(x,y) + 2*(1 + A)*(-2 + wx(x,y))*(-1 + wx(x,y) + wy(x,y))* dxpy(x,y)))))/(6*(1 + A)*pi* wx(x,y)*(-2 + wx(x,y) + wy(x,y))*(-1 + wx(x,y) + wy(x,y))^(3/2)* T(x,y)) ) ;
      end
    end
  end

















    phydro();
    %pause(0.033);

  function r = phydro()
    frame_active = (mod(t-1,int16(1/DT/SECONDS_PER_LIGHTYEAR/30)) == 0);
    if VIDEO & ~frame_active 
      return
    end
    if VIDEO
      disp(t);
    end
    plotsx = 4;
    plotsy = 2;

    subplot(plotsy, plotsx, 1);
    image((T-TEMP_OFFSET) * 256 + 128); 
    
    % if (VIDEO) & (frame_active)
    %   frame = getframe;
    %   writeVideo(vt,frame);
    % end

    subplot(plotsy, plotsx, 2);
    image(vx*128 + 128);
    % if (VIDEO) & (frame_active)
    %   frame = getframe;
    %   writeVideo(vvx,frame);
    % end

    subplot(plotsy, plotsx, 3);
    image(vy*128 + 128);
    % if (VIDEO) & (frame_active)
    %   frame = getframe;
    %   writeVideo(vvy,frame);
    % end

    subplot(plotsy, plotsx, 4);
    image(sqrt(vx.*vx + vy.*vy)*128 + 128);
    % if (VIDEO) & (frame_active)
    %   frame = getframe;
    %   writeVideo(vv,frame);
    % end

    subplot(plotsy, plotsx, 5);
    image(px*8 + 128);

    subplot(plotsy, plotsx, 6);
    image(py*8 + 128);
    truesize(fig);
    if (VIDEO) & (frame_active)
      frame = getframe(gcf);
      writeVideo(vt,frame);
    end
    
    % subplot(plotsy, plotsx, 9);
    % image(dtdtvx*8 + 128);
    % subplot(plotsy, plotsx, 10);
    % image(dtdtvy*8 + 128);
    % subplot(plotsy, plotsx, 11);
    % image(Pi*2 + 128);
    % subplot(plotsy, plotsx, 12);
    % image(B*2 + 128);
    return;

    subplot(2, 6, 2);
    image(-Pi*4); 
    subplot(2, 6, 3);
    image(B*4); 
    subplot(2, 6, 4);
    image(-B*4); 
    subplot(2, 6, 5);
    image(px); 
    subplot(2, 6, 6);
    image(py); 
    subplot(2, 6, 9);
    image(dtdtvx*10000); 
    subplot(2, 6, 10);
    image(dtdtvy*100000); 
    subplot(2, 6, 11);
    image(dtpx); 
    subplot(2, 6, 12);
    image(dtpy); 
    
    %subplot(1, 6, 6);
    %image(T*256); 
    truesize(fig);
  end

  % [-1 0 1
  function r = dx2(arg)
    r = zeros(GRID_SIZE, GRID_SIZE);
    for y = 1:GRID_SIZE
      r(1,y) = - arg(GRID_SIZE,y)/2 + arg(2, y)/2;
    end
    for y = 1:GRID_SIZE
      r(GRID_SIZE,y) =  -arg(GRID_SIZE-1,y)/2 + arg(1, y)/2;
    end
    for x = 2:GRID_SIZE-1
      for y = 1:GRID_SIZE
        r(x, y) = (arg(x+1,y) - arg(x-1,y))/2;   
      end
    end 
    r = r / DX;
  end
  function r = dy2(arg)
    r = zeros(GRID_SIZE, GRID_SIZE);
    for x = 1:GRID_SIZE
      r(x,1) = - arg(x,GRID_SIZE)/2 + arg(x,2)/2;
    end
    for x = 1:GRID_SIZE
      r(x,GRID_SIZE) =  -arg(x,GRID_SIZE-1)/2 + arg(x,1)/2;
    end
    for x = 1:GRID_SIZE
      for y = 2:GRID_SIZE-1
        r(x, y) = (arg(x,y+1) - arg(x,y-1))/2;   
      end
    end 
    r = r / DX;
  end
  function r = dx(arg)
    r = d1 * arg;
    r = r / DX;
  end
  function r = dy(arg)
    r = arg * (d1.');
    r = r / DX;
  end
  function r = dxdx(arg)
    r = d2 * arg;
    r = r / DX / DX;
  end
  function r = dydy(arg)
    r = arg * (d2.');
    r = r / DX / DX;
  end
  function r = dxdx2(arg)
    r = zeros(GRID_SIZE, GRID_SIZE);
    for y = 1:GRID_SIZE
      r(1,y) = arg(GRID_SIZE,y) + arg(2, y) - 2*arg(1,y);
    end
    for y = 1:GRID_SIZE
      r(GRID_SIZE,y) =  arg(GRID_SIZE-1,y) + arg(1, y) - 2*arg(GRID_SIZE,y);
    end
    for x = 2:GRID_SIZE-1
      for y = 1:GRID_SIZE
        r(x, y) = arg(x+1,y) + arg(x-1,y) - 2*arg(x,y);
      end
    end 
    r = r / DX / DX;
  end
  function r = dydy2(arg)
    r = zeros(GRID_SIZE, GRID_SIZE);
    for x = 1:GRID_SIZE
      r(x,1) = arg(x,GRID_SIZE) + arg(x,2) - 2*arg(x,1);
    end
    for x = 1:GRID_SIZE
      r(x,GRID_SIZE) =  arg(x,GRID_SIZE-1) + arg(x,1) - 2*arg(x,GRID_SIZE);
    end
    for x = 1:GRID_SIZE
      for y = 2:GRID_SIZE-1
        r(x, y) = arg(x,y+1) + arg(x,y-1) - 2*arg(x,y);   
      end
    end 
    r = r / DX / DX;
  end

end


